HerbPt1 16S Figures

require(phyloseq)

Loading required package: phyloseq

require(tidyverse)

Loading required package: tidyverse
Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.1 ──
✓ ggplot2 3.3.5     ✓ purrr   0.3.4
✓ tibble  3.1.6     ✓ dplyr   1.0.8
✓ tidyr   1.2.0     ✓ stringr 1.4.0
✓ readr   2.1.2     ✓ forcats 0.5.1
Warning: package ‘tidyr’ was built under R version 4.1.2
Warning: package ‘readr’ was built under R version 4.1.2
Warning: package ‘dplyr’ was built under R version 4.1.2
── Conflicts ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

require(reshape2)

Loading required package: reshape2

Attaching package: ‘reshape2’

The following object is masked from ‘package:tidyr’:

    smiths

require(dplyr)
require(ggplot2)
require(vegan)

Loading required package: vegan
Warning: package ‘vegan’ was built under R version 4.1.2
Loading required package: permute
Warning: package ‘permute’ was built under R version 4.1.2
Loading required package: lattice
This is vegan 2.6-2

Load data order, factors, and create a mode (chemical, hand, non-treated) column.

ps_dmn <- readRDS("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/data/PhyloseqObjects/16S/DMN_ests_16S.Rdata")
sample_data(ps_dmn)$Herbicide <- factor(sample_data(ps_dmn)$Herbicide, levels = c("Aatrex", "Clarity", "Hand","Non-Treated","Roundup Powermax"))
sample_data(ps_dmn)$herb_time<-paste(sample_data(ps_dmn)$Herbicide, sample_data(ps_dmn)$Time, sep = "_")

#regroup all chemical treatments together and rerun betadiv calcs within group. 
sample_data(ps_dmn)$Mode<-sample_data(ps_dmn)$Herbicide

index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Chemical", "Chemical", "Chemical", "Hand", "Non-Treated")

sample_data(ps_dmn)$Mode<- as.factor(values[match(sample_data(ps_dmn)$Mode, index)])


index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Dicamba", "Glyphosate", "Atrazine-Mesotrione", "Handweeded", "Non-Treated")

sample_data(ps_dmn)$Herbicide <- as.factor(values[match(sample_data(ps_dmn)$Herbicide, index)])



ps_rare <- readRDS("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/data/PhyloseqObjects/16S/HerbPt1_rare_16S.Rdata")
sample_data(ps_rare)$Herbicide <- factor(sample_data(ps_rare)$Herbicide, levels = c("Aatrex", "Clarity", "Hand","Non-Treated","Roundup Powermax"))
sample_data(ps_rare)$herb_time<-paste(sample_data(ps_rare)$Herbicide, sample_data(ps_rare)$Time, sep = "_")


#regroup all chemical treatments together and rerun betadiv calcs within group. 
sample_data(ps_rare)$Mode<-sample_data(ps_rare)$Herbicide

index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Chemical", "Chemical", "Chemical", "Hand", "Non-Treated")

sample_data(ps_rare)$Mode<- as.factor(values[match(sample_data(ps_rare)$Mode, index)])

index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Dicamba", "Glyphosate", "Atrazine-Mesotrione", "Handweeded", "Non-Treated")

sample_data(ps_rare)$Herbicide <- as.factor(values[match(sample_data(ps_rare)$Herbicide, index)])

ps_trans <- readRDS("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/data/PhyloseqObjects/16S/HerbPt1_hel_trans_16S.Rdata")
sample_data(ps_trans)$Herbicide <- factor(sample_data(ps_trans)$Herbicide, levels = c("Aatrex", "Clarity", "Hand","Non-Treated","Roundup Powermax"))
sample_data(ps_trans)$herb_time<-paste(sample_data(ps_trans)$Herbicide, sample_data(ps_trans)$Time, sep = "_")

#regroup all chemical treatments together and rerun betadiv calcs within group. 
sample_data(ps_trans)$Mode<-sample_data(ps_trans)$Herbicide

index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Chemical", "Chemical", "Chemical", "Hand", "Non-Treated")

sample_data(ps_trans)$Mode<- as.factor(values[match(sample_data(ps_trans)$Mode, index)])

index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Dicamba", "Glyphosate", "Atrazine-Mesotrione", "Handweeded", "Non-Treated")

sample_data(ps_trans)$Herbicide <- as.factor(values[match(sample_data(ps_trans)$Herbicide, index)])

create alpha diversity tables

alpha_div_md$Herbicide

  [1] Dicamba             Dicamba             Dicamba             Handweeded          Handweeded          Handweeded          Non-Treated         Non-Treated        
  [9] Non-Treated         Atrazine-Mesotrione Atrazine-Mesotrione Glyphosate          Glyphosate          Glyphosate          Glyphosate          Glyphosate         
 [17] Atrazine-Mesotrione Atrazine-Mesotrione Atrazine-Mesotrione Handweeded          Handweeded          Dicamba             Dicamba             Dicamba            
 [25] Non-Treated         Non-Treated         Non-Treated         Non-Treated         Non-Treated         Non-Treated         Glyphosate          Glyphosate         
 [33] Glyphosate          Dicamba             Dicamba             Dicamba             Atrazine-Mesotrione Atrazine-Mesotrione Atrazine-Mesotrione Handweeded         
 [41] Handweeded          Handweeded          Glyphosate          Glyphosate          Glyphosate          Handweeded          Handweeded          Non-Treated        
 [49] Non-Treated         Non-Treated         Dicamba             Dicamba             Dicamba             Atrazine-Mesotrione Atrazine-Mesotrione Atrazine-Mesotrione
 [57] Dicamba             Dicamba             Handweeded          Handweeded          Handweeded          Non-Treated         Non-Treated         Non-Treated        
 [65] Atrazine-Mesotrione Atrazine-Mesotrione Atrazine-Mesotrione Glyphosate          Glyphosate          Glyphosate          Glyphosate          Glyphosate         
 [73] Glyphosate          Atrazine-Mesotrione Atrazine-Mesotrione Atrazine-Mesotrione Handweeded          Handweeded          Dicamba             Dicamba            
 [81] Dicamba             Non-Treated         Non-Treated         Non-Treated         Non-Treated         Non-Treated         Glyphosate          Glyphosate         
 [89] Dicamba             Dicamba             Dicamba             Atrazine-Mesotrione Atrazine-Mesotrione Atrazine-Mesotrione Handweeded          Handweeded         
 [97] Handweeded          Glyphosate          Glyphosate          Handweeded          Handweeded          Handweeded          Non-Treated         Non-Treated        
[105] Non-Treated         Dicamba             Dicamba             Dicamba             Atrazine-Mesotrione Atrazine-Mesotrione Dicamba             Dicamba            
[113] Dicamba             Handweeded          Handweeded          Handweeded          Non-Treated         Non-Treated         Atrazine-Mesotrione Atrazine-Mesotrione
[121] Atrazine-Mesotrione Glyphosate          Glyphosate          Glyphosate          Glyphosate          Glyphosate          Atrazine-Mesotrione Atrazine-Mesotrione
[129] Atrazine-Mesotrione Handweeded          Handweeded          Handweeded          Dicamba             Dicamba             Dicamba             Non-Treated        
[137] Non-Treated         Non-Treated         Non-Treated         Non-Treated         Glyphosate          Glyphosate          Glyphosate          Dicamba            
[145] Dicamba             Dicamba             Atrazine-Mesotrione Atrazine-Mesotrione Atrazine-Mesotrione Handweeded          Handweeded          Glyphosate         
[153] Glyphosate          Glyphosate          Handweeded          Handweeded          Handweeded          Non-Treated         Non-Treated         Non-Treated        
[161] Dicamba             Dicamba             Atrazine-Mesotrione Atrazine-Mesotrione Atrazine-Mesotrione
Levels: Non-Treated Handweeded Atrazine-Mesotrione Dicamba Glyphosate

Shannon Div plots - no significant differences among herbicide treatments at any of the three time points

ggplot(data = alpha_div_md, aes(Herbicide, Shannon, color= Herbicide)) + facet_grid(. ~ Time) + geom_boxplot() + theme_classic() + theme(axis.text.x = element_text(angle = 45, hjust = 1) )

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_Shannon.pdf")

Saving 7.29 x 4.51 in image

aov_t1<-aov(Chao1 ~ Herbicide, data = alpha_div_md[alpha_div_md$Time == "T1",])
plot(aov_t1$residuals)

summary(aov_t1)

            Df  Sum Sq Mean Sq F value Pr(>F)
Herbicide    4   66537   16634   0.099  0.982
Residuals   51 8553352  167713               

aov_t2<-aov(Chao1 ~ Herbicide, data = alpha_div_md[alpha_div_md$Time == "T2",])
plot(aov_t2$residuals)

summary(aov_t2)

            Df   Sum Sq Mean Sq F value Pr(>F)
Herbicide    4   652901  163225   0.779  0.545
Residuals   49 10272922  209651               

aov_t3<-aov(Chao1 ~ Herbicide, data = alpha_div_md[alpha_div_md$Time == "T3",])
plot(aov_t3$residuals)

summary(aov_t3)

            Df   Sum Sq Mean Sq F value Pr(>F)
Herbicide    4   433435  108359   0.465  0.761
Residuals   50 11641251  232825               

remove outliers and rare reads with less than 2 total reads

ps_dmn <-  subset_samples(ps_dmn, sample_names(ps_rare) != "G166SG")

ps_rare <-  subset_samples(ps_rare, sample_names(ps_rare) != "G166SG")
ps_rare_sub<-prune_taxa(taxa_sums(ps_rare) > 2, ps_rare)

ps_trans_sub<-prune_taxa(taxa_sums(ps_trans) > 0.05, ps_trans)

ordinations and adonis testing with three separate objects (i.e., dmn, rarefied, transformed). Rare taxa are removed from rarefied and transfomred to sucessfully ordinate. At this point, the transformed data will not ordinate.

ord_dmn<-ordinate(physeq = ps_dmn, method = "NMDS", distance = "bray", k=3, trymax= 300, maxit=1000)

Run 0 stress 0.1117941 
Run 1 stress 0.1115908 
... New best solution
... Procrustes: rmse 0.01112287  max resid 0.137285 
Run 2 stress 0.1118138 
... Procrustes: rmse 0.01224305  max resid 0.1486074 
Run 3 stress 0.1117997 
... Procrustes: rmse 0.01238183  max resid 0.1520467 
Run 4 stress 0.1118133 
... Procrustes: rmse 0.01165924  max resid 0.1445305 
Run 5 stress 0.1117717 
... Procrustes: rmse 0.004992797  max resid 0.05385656 
Run 6 stress 0.1115916 
... Procrustes: rmse 0.0003952995  max resid 0.003519588 
... Similar to previous best
Run 7 stress 0.1117649 
... Procrustes: rmse 0.004519199  max resid 0.05419881 
Run 8 stress 0.1115917 
... Procrustes: rmse 0.001064847  max resid 0.01208441 
Run 9 stress 0.1117663 
... Procrustes: rmse 0.00470262  max resid 0.0554819 
Run 10 stress 0.1159638 
Run 11 stress 0.1117987 
... Procrustes: rmse 0.01229483  max resid 0.1505899 
Run 12 stress 0.1118121 
... Procrustes: rmse 0.01145765  max resid 0.1419449 
Run 13 stress 0.1118606 
... Procrustes: rmse 0.008173993  max resid 0.09623366 
Run 14 stress 0.111798 
... Procrustes: rmse 0.01121996  max resid 0.14002 
Run 15 stress 0.1118121 
... Procrustes: rmse 0.01340603  max resid 0.1643183 
Run 16 stress 0.1115947 
... Procrustes: rmse 0.001573293  max resid 0.01899481 
Run 17 stress 0.1118639 
... Procrustes: rmse 0.008635792  max resid 0.1014058 
Run 18 stress 0.1115911 
... Procrustes: rmse 0.0003000347  max resid 0.002739309 
... Similar to previous best
Run 19 stress 0.1115954 
... Procrustes: rmse 0.001243378  max resid 0.01498226 
Run 20 stress 0.1117972 
... Procrustes: rmse 0.01216683  max resid 0.1493349 
*** Solution reached

ord_rare<-ordinate(physeq = ps_rare_sub, method = "NMDS", distance = "bray", k=3, trymax= 300, maxit=1000)

Square root transformation
Wisconsin double standardization
Run 0 stress 0.1959743 
Run 1 stress 0.1955267 
... New best solution
... Procrustes: rmse 0.01800803  max resid 0.1123914 
Run 2 stress 0.1955394 
... Procrustes: rmse 0.007006262  max resid 0.05765258 
Run 3 stress 0.1958487 
... Procrustes: rmse 0.01579267  max resid 0.1056545 
Run 4 stress 0.1955945 
... Procrustes: rmse 0.01295785  max resid 0.1494722 
Run 5 stress 0.1955967 
... Procrustes: rmse 0.01293347  max resid 0.1501704 
Run 6 stress 0.1969796 
Run 7 stress 0.1955973 
... Procrustes: rmse 0.01337124  max resid 0.1534108 
Run 8 stress 0.1955967 
... Procrustes: rmse 0.01339074  max resid 0.1533846 
Run 9 stress 0.1959732 
... Procrustes: rmse 0.01729624  max resid 0.1080525 
Run 10 stress 0.1956108 
... Procrustes: rmse 0.01142375  max resid 0.1048192 
Run 11 stress 0.1956435 
... Procrustes: rmse 0.01260475  max resid 0.1537692 
Run 12 stress 0.1955174 
... New best solution
... Procrustes: rmse 0.008122313  max resid 0.05712863 
Run 13 stress 0.1958337 
... Procrustes: rmse 0.01090413  max resid 0.1036971 
Run 14 stress 0.1967865 
Run 15 stress 0.1955939 
... Procrustes: rmse 0.01324124  max resid 0.1490696 
Run 16 stress 0.196786 
Run 17 stress 0.1955234 
... Procrustes: rmse 0.007498469  max resid 0.05695487 
Run 18 stress 0.195582 
... Procrustes: rmse 0.01173812  max resid 0.07111958 
Run 19 stress 0.1958207 
... Procrustes: rmse 0.009376553  max resid 0.1035924 
Run 20 stress 0.1967535 
Run 21 stress 0.1956211 
... Procrustes: rmse 0.01156007  max resid 0.1369555 
Run 22 stress 0.195582 
... Procrustes: rmse 0.01170171  max resid 0.07088219 
Run 23 stress 0.1955424 
... Procrustes: rmse 0.009362655  max resid 0.06422712 
Run 24 stress 0.1954828 
... New best solution
... Procrustes: rmse 0.006018466  max resid 0.04375112 
Run 25 stress 0.1955168 
... Procrustes: rmse 0.00608034  max resid 0.04463714 
Run 26 stress 0.1965865 
Run 27 stress 0.196778 
Run 28 stress 0.1955925 
... Procrustes: rmse 0.01221595  max resid 0.1482877 
Run 29 stress 0.1955944 
... Procrustes: rmse 0.01222552  max resid 0.1482843 
Run 30 stress 0.1958206 
... Procrustes: rmse 0.01151099  max resid 0.1048252 
Run 31 stress 0.1956224 
... Procrustes: rmse 0.01254164  max resid 0.134756 
Run 32 stress 0.1971977 
Run 33 stress 0.1969766 
Run 34 stress 0.195833 
... Procrustes: rmse 0.01382711  max resid 0.1051923 
Run 35 stress 0.1954885 
... Procrustes: rmse 0.001224376  max resid 0.01288981 
Run 36 stress 0.1954874 
... Procrustes: rmse 0.001327997  max resid 0.01348466 
Run 37 stress 0.1967515 
Run 38 stress 0.1958228 
... Procrustes: rmse 0.01137911  max resid 0.1048475 
Run 39 stress 0.197019 
Run 40 stress 0.1958336 
... Procrustes: rmse 0.01400106  max resid 0.1052911 
Run 41 stress 0.1955898 
... Procrustes: rmse 0.01087936  max resid 0.1192967 
Run 42 stress 0.1954861 
... Procrustes: rmse 0.001352308  max resid 0.01107595 
Run 43 stress 0.1984585 
Run 44 stress 0.1955182 
... Procrustes: rmse 0.005903524  max resid 0.04272373 
Run 45 stress 0.1970057 
Run 46 stress 0.1965857 
Run 47 stress 0.1974123 
Run 48 stress 0.1956422 
... Procrustes: rmse 0.01285919  max resid 0.1537516 
Run 49 stress 0.1956489 
... Procrustes: rmse 0.01316701  max resid 0.1561775 
Run 50 stress 0.195849 
... Procrustes: rmse 0.01454653  max resid 0.1054922 
Run 51 stress 0.195821 
... Procrustes: rmse 0.01228135  max resid 0.105193 
Run 52 stress 0.1978668 
Run 53 stress 0.1958338 
... Procrustes: rmse 0.01282964  max resid 0.10528 
Run 54 stress 0.1955268 
... Procrustes: rmse 0.003679333  max resid 0.03000048 
Run 55 stress 0.1971375 
Run 56 stress 0.1956175 
... Procrustes: rmse 0.01127198  max resid 0.1151186 
Run 57 stress 0.1956146 
... Procrustes: rmse 0.01061111  max resid 0.1067233 
Run 58 stress 0.1958342 
... Procrustes: rmse 0.01302074  max resid 0.1053564 
Run 59 stress 0.1972164 
Run 60 stress 0.196753 
Run 61 stress 0.197357 
Run 62 stress 0.1956426 
... Procrustes: rmse 0.01282273  max resid 0.1539347 
Run 63 stress 0.1956151 
... Procrustes: rmse 0.01074056  max resid 0.1084906 
Run 64 stress 0.195649 
... Procrustes: rmse 0.01317392  max resid 0.1561167 
Run 65 stress 0.1959754 
... Procrustes: rmse 0.01604176  max resid 0.1104391 
Run 66 stress 0.1954887 
... Procrustes: rmse 0.001712985  max resid 0.01764067 
Run 67 stress 0.1969292 
Run 68 stress 0.1967917 
Run 69 stress 0.1955019 
... Procrustes: rmse 0.003005991  max resid 0.01927756 
Run 70 stress 0.196758 
Run 71 stress 0.1972331 
Run 72 stress 0.1971941 
Run 73 stress 0.1954839 
... Procrustes: rmse 0.000538266  max resid 0.005283253 
... Similar to previous best
*** Solution reached

full_ord_rare<-ggordiplots::gg_ordiplot(ord = ord_rare, groups = data.frame(sample_data(ps_rare))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)

full_ord_rare$plot + theme_classic()

ord_transformed<-ordinate(physeq = ps_trans_sub, method = "NMDS", distance = "bray", trymax= 300, maxit=1000)

Run 0 stress 0.07436383 
Run 1 stress 0.07411606 
... New best solution
... Procrustes: rmse 0.03117035  max resid 0.2498244 
Run 2 stress 0.07773971 
Run 3 stress 0.07217597 
... New best solution
... Procrustes: rmse 0.05313851  max resid 0.4765044 
Run 4 stress 0.07455535 
Run 5 stress 0.07413861 
Run 6 stress 0.07821068 
Run 7 stress 0.07475795 
Run 8 stress 0.07634371 
Run 9 stress 0.07390696 
Run 10 stress 0.07475868 
Run 11 stress 0.07371342 
Run 12 stress 0.07653065 
Run 13 stress 0.07356773 
Run 14 stress 0.07214992 
... New best solution
... Procrustes: rmse 0.001443487  max resid 0.01837558 
Run 15 stress 0.07478238 
Run 16 stress 0.07401673 
Run 17 stress 0.07564706 
Run 18 stress 0.07376447 
Run 19 stress 0.0763201 
Run 20 stress 0.08031617 
Run 21 stress 0.07350122 
Run 22 stress 0.07472829 
Run 23 stress 0.07455301 
Run 24 stress 0.07224441 
... Procrustes: rmse 0.002650218  max resid 0.02493022 
Run 25 stress 0.07343235 
Run 26 stress 0.08340946 
Run 27 stress 0.07371382 
Run 28 stress 0.07809551 
Run 29 stress 0.07565057 
Run 30 stress 0.07604995 
Run 31 stress 0.07400801 
Run 32 stress 0.07455409 
Run 33 stress 0.08407341 
Run 34 stress 0.08405253 
Run 35 stress 0.08494465 
Run 36 stress 0.07207254 
... New best solution
... Procrustes: rmse 0.02923975  max resid 0.2677599 
Run 37 stress 0.07433151 
Run 38 stress 0.08031236 
Run 39 stress 0.07401128 
Run 40 stress 0.07630296 
Run 41 stress 0.07504806 
Run 42 stress 0.07353092 
Run 43 stress 0.07210215 
... Procrustes: rmse 0.001996828  max resid 0.0169429 
Run 44 stress 0.07290392 
Run 45 stress 0.07456734 
Run 46 stress 0.4158745 
Run 47 stress 0.07468176 
Run 48 stress 0.07372598 
Run 49 stress 0.07314001 
Run 50 stress 0.07431572 
Run 51 stress 0.07217868 
... Procrustes: rmse 0.02929952  max resid 0.2671141 
Run 52 stress 0.07401337 
Run 53 stress 0.07419621 
Run 54 stress 0.0761684 
Run 55 stress 0.07390179 
Run 56 stress 0.08495755 
Run 57 stress 0.07457704 
Run 58 stress 0.07372457 
Run 59 stress 0.07411643 
Run 60 stress 0.07810726 
Run 61 stress 0.07435209 
Run 62 stress 0.07565289 
Run 63 stress 0.07399032 
Run 64 stress 0.08507366 
Run 65 stress 0.07419598 
Run 66 stress 0.08221774 
Run 67 stress 0.07388912 
Run 68 stress 0.08249078 
Run 69 stress 0.07401304 
Run 70 stress 0.07388093 
Run 71 stress 0.07399428 
Run 72 stress 0.08033213 
Run 73 stress 0.07371329 
Run 74 stress 0.07379163 
Run 75 stress 0.08424884 
Run 76 stress 0.07382315 
Run 77 stress 0.07495688 
Run 78 stress 0.07370127 
Run 79 stress 0.07475771 
Run 80 stress 0.07215002 
... Procrustes: rmse 0.02929661  max resid 0.2675789 
Run 81 stress 0.07453389 
Run 82 stress 0.07291297 
Run 83 stress 0.08293632 
Run 84 stress 0.07720927 
Run 85 stress 0.07438212 
Run 86 stress 0.07390787 
Run 87 stress 0.07434478 
Run 88 stress 0.07433927 
Run 89 stress 0.08539586 
Run 90 stress 0.08382492 
Run 91 stress 0.07398725 
Run 92 stress 0.07407376 
Run 93 stress 0.07353074 
Run 94 stress 0.07432943 
Run 95 stress 0.07207262 
... Procrustes: rmse 0.0004636339  max resid 0.004175306 
... Similar to previous best
*** Solution reached

Adonis testing of herbicide treatments by time point

ps_adonis<-function(physeq){
  otu_tab<-data.frame(phyloseq::otu_table(physeq))
  md_tab<-data.frame(phyloseq::sample_data(physeq))
    if(taxa_are_rows(physeq)== T){
       physeq_dist<-parallelDist::parDist(as.matrix(t(otu_tab)), method = "bray")}
            else{physeq_dist<-parallelDist::parDist(as.matrix(otu_tab), method = "bray")}
  print(anova(vegan::betadisper(physeq_dist, md_tab$Herbicide)))
  vegan::adonis(physeq_dist ~ Herbicide * Time + Total_Weed_Veg , data = md_tab, permutations = 1000)
}

remove one sample with no vegetation measurement.

ps_rare_sub_57<-subset_samples(ps_rare_sub, sample_names(ps_rare_sub) != "G065SG")
ps_adonis(ps_rare_sub_57)

ps_dmn_57<-subset_samples(ps_dmn, sample_names(ps_dmn) != "G065SG")
ps_adonis(ps_dmn_57)

Ordination plots DMN

ord_t1_dmn<-ordinate(physeq = subset_samples(ps_dmn, Time=="T1"), method = "NMDS", distance = "bray", k=3, trymax= 100)

Run 0 stress 0.0860376 
Run 1 stress 0.08603806 
... Procrustes: rmse 0.0001239324  max resid 0.0007295483 
... Similar to previous best
Run 2 stress 0.0983852 
Run 3 stress 0.08674017 
Run 4 stress 0.09974458 
Run 5 stress 0.08662921 
Run 6 stress 0.08662975 
Run 7 stress 0.08674026 
Run 8 stress 0.08603767 
... Procrustes: rmse 2.019457e-05  max resid 0.0001271086 
... Similar to previous best
Run 9 stress 0.1024486 
Run 10 stress 0.0866297 
Run 11 stress 0.0866588 
Run 12 stress 0.0860374 
... New best solution
... Procrustes: rmse 0.0008541763  max resid 0.004577351 
... Similar to previous best
Run 13 stress 0.08674053 
Run 14 stress 0.08662909 
Run 15 stress 0.0997284 
Run 16 stress 0.08603683 
... New best solution
... Procrustes: rmse 0.0002639494  max resid 0.001218145 
... Similar to previous best
Run 17 stress 0.08674004 
Run 18 stress 0.08603779 
... Procrustes: rmse 0.0006511547  max resid 0.003685441 
... Similar to previous best
Run 19 stress 0.08673995 
Run 20 stress 0.08673991 
*** Solution reached

T1_dmn<-ggordiplots::gg_ordiplot(ord = ord_t1_dmn, groups = data.frame(sample_data(subset_samples(ps_dmn, Time == "T1")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)

T1_dmn$plot + theme_classic() + xlim(-0.4, 0.4) + ylim(-0.4, 0.4)

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_dmn_T1.pdf")

Saving 7.29 x 4.51 in image

ord_t2_dmn<-ordinate(physeq = subset_samples(ps_dmn, Time=="T2"), method = "NMDS", distance = "bray", k=3, trymax= 100)

Run 0 stress 0.1104159 
Run 1 stress 0.1104153 
... New best solution
... Procrustes: rmse 0.0003571464  max resid 0.002189988 
... Similar to previous best
Run 2 stress 0.1095397 
... New best solution
... Procrustes: rmse 0.03309861  max resid 0.2348476 
Run 3 stress 0.1096257 
... Procrustes: rmse 0.005714701  max resid 0.03065104 
Run 4 stress 0.1095261 
... New best solution
... Procrustes: rmse 0.0006951312  max resid 0.003638642 
... Similar to previous best
Run 5 stress 0.1104161 
Run 6 stress 0.1122104 
Run 7 stress 0.1095262 
... Procrustes: rmse 0.0001081361  max resid 0.0004776 
... Similar to previous best
Run 8 stress 0.1122142 
Run 9 stress 0.1095263 
... Procrustes: rmse 0.0006652566  max resid 0.003388802 
... Similar to previous best
Run 10 stress 0.1095263 
... Procrustes: rmse 0.0005610351  max resid 0.003406831 
... Similar to previous best
Run 11 stress 0.1108789 
Run 12 stress 0.1095267 
... Procrustes: rmse 0.0008601728  max resid 0.004694253 
... Similar to previous best
Run 13 stress 0.1104149 
Run 14 stress 0.1122102 
Run 15 stress 0.1104149 
Run 16 stress 0.1108808 
Run 17 stress 0.1104147 
Run 18 stress 0.1105647 
Run 19 stress 0.1108803 
Run 20 stress 0.1108777 
*** Solution reached

T2_dmn<-ggordiplots::gg_ordiplot(ord = ord_t2_dmn, groups = data.frame(sample_data(subset_samples(ps_dmn, Time == "T2")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1) 

T2_dmn$plot + theme_classic()+ xlim(-0.4, 0.4) + ylim(-0.4, 0.4)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_dmn_T2.pdf")

Saving 7.29 x 4.51 in image

ord_t3_dmn<-ordinate(physeq = subset_samples(ps_dmn, Time=="T3"), method = "NMDS", distance = "bray", k=3, trymax= 100)

Run 0 stress 0.09546215 
Run 1 stress 0.09706147 
Run 2 stress 0.09546159 
... New best solution
... Procrustes: rmse 0.001277492  max resid 0.007259895 
... Similar to previous best
Run 3 stress 0.0955498 
... Procrustes: rmse 0.03810422  max resid 0.2424179 
Run 4 stress 0.09550233 
... Procrustes: rmse 0.03778166  max resid 0.2444873 
Run 5 stress 0.09546332 
... Procrustes: rmse 0.001100842  max resid 0.005583535 
... Similar to previous best
Run 6 stress 0.09732995 
Run 7 stress 0.09546077 
... New best solution
... Procrustes: rmse 0.0007947735  max resid 0.00456194 
... Similar to previous best
Run 8 stress 0.09546428 
... Procrustes: rmse 0.001848437  max resid 0.009343808 
... Similar to previous best
Run 9 stress 0.09712155 
Run 10 stress 0.09549825 
... Procrustes: rmse 0.03714564  max resid 0.2419844 
Run 11 stress 0.09549948 
... Procrustes: rmse 0.0373624  max resid 0.2433029 
Run 12 stress 0.09546252 
... Procrustes: rmse 0.0005053034  max resid 0.002683928 
... Similar to previous best
Run 13 stress 0.09569417 
... Procrustes: rmse 0.03993933  max resid 0.2524409 
Run 14 stress 0.09546543 
... Procrustes: rmse 0.00135589  max resid 0.007233676 
... Similar to previous best
Run 15 stress 0.0971071 
Run 16 stress 0.09549812 
... Procrustes: rmse 0.03722254  max resid 0.2427672 
Run 17 stress 0.0970629 
Run 18 stress 0.09546061 
... New best solution
... Procrustes: rmse 0.000334493  max resid 0.00147327 
... Similar to previous best
Run 19 stress 0.09570193 
... Procrustes: rmse 0.04006183  max resid 0.2508731 
Run 20 stress 0.09549861 
... Procrustes: rmse 0.03739865  max resid 0.2432734 
*** Solution reached

T3_dmn<-ggordiplots::gg_ordiplot(ord = ord_t3_dmn, groups = data.frame(sample_data(subset_samples(ps_dmn, Time == "T3")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1) 

T3_dmn$plot + theme_classic()+ xlim(-0.4, 0.4) + ylim(-0.4, 0.4)

Warning: Removed 1 rows containing missing values (geom_point).

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_dmn_T3.pdf")

Saving 7.29 x 4.51 in image
Warning: Removed 1 rows containing missing values (geom_point).

Ordination plots rarefied

ord_t1_rare<-ordinate(physeq = subset_samples(ps_rare, Time=="T1"), method = "NMDS", distance = "bray", k=3, trymax= 100)

Square root transformation
Wisconsin double standardization
Run 0 stress 0.1810491 
Run 1 stress 0.1812643 
... Procrustes: rmse 0.02583058  max resid 0.1092403 
Run 2 stress 0.1814142 
... Procrustes: rmse 0.04648293  max resid 0.1540837 
Run 3 stress 0.1824085 
Run 4 stress 0.1854315 
Run 5 stress 0.1817696 
Run 6 stress 0.1813795 
... Procrustes: rmse 0.01903201  max resid 0.1229428 
Run 7 stress 0.1873045 
Run 8 stress 0.1810319 
... New best solution
... Procrustes: rmse 0.03213373  max resid 0.1238225 
Run 9 stress 0.1815254 
... Procrustes: rmse 0.07330822  max resid 0.260521 
Run 10 stress 0.181525 
... Procrustes: rmse 0.07331803  max resid 0.2605261 
Run 11 stress 0.1815251 
... Procrustes: rmse 0.07331039  max resid 0.2605071 
Run 12 stress 0.1812636 
... Procrustes: rmse 0.02331833  max resid 0.1081706 
Run 13 stress 0.1873031 
Run 14 stress 0.1809659 
... New best solution
... Procrustes: rmse 0.04119361  max resid 0.1474758 
Run 15 stress 0.1822765 
Run 16 stress 0.1864799 
Run 17 stress 0.1815253 
Run 18 stress 0.1837641 
Run 19 stress 0.1814576 
... Procrustes: rmse 0.05589707  max resid 0.2469422 
Run 20 stress 0.1836188 
Run 21 stress 0.1819126 
Run 22 stress 0.1821857 
Run 23 stress 0.1810495 
... Procrustes: rmse 0.01875174  max resid 0.1152783 
Run 24 stress 0.1822761 
Run 25 stress 0.1808423 
... New best solution
... Procrustes: rmse 0.03509173  max resid 0.1424772 
Run 26 stress 0.1815222 
Run 27 stress 0.1814211 
Run 28 stress 0.1881278 
Run 29 stress 0.1831937 
Run 30 stress 0.1815254 
Run 31 stress 0.1815254 
Run 32 stress 0.1877224 
Run 33 stress 0.1823821 
Run 34 stress 0.1816035 
Run 35 stress 0.1828841 
Run 36 stress 0.1811594 
... Procrustes: rmse 0.02182468  max resid 0.1078205 
Run 37 stress 0.181438 
Run 38 stress 0.182276 
Run 39 stress 0.1824766 
Run 40 stress 0.1875966 
Run 41 stress 0.182182 
Run 42 stress 0.1878771 
Run 43 stress 0.1872815 
Run 44 stress 0.1885431 
Run 45 stress 0.1821561 
Run 46 stress 0.1815255 
Run 47 stress 0.1816037 
Run 48 stress 0.1810502 
... Procrustes: rmse 0.02488083  max resid 0.1084949 
Run 49 stress 0.1842028 
Run 50 stress 0.1831935 
Run 51 stress 0.181416 
Run 52 stress 0.1825128 
Run 53 stress 0.1821252 
Run 54 stress 0.1882 
Run 55 stress 0.1858117 
Run 56 stress 0.1910572 
Run 57 stress 0.1814036 
Run 58 stress 0.1815253 
Run 59 stress 0.1837646 
Run 60 stress 0.1855356 
Run 61 stress 0.1886207 
Run 62 stress 0.1811714 
... Procrustes: rmse 0.02206901  max resid 0.1101655 
Run 63 stress 0.1813283 
... Procrustes: rmse 0.02669062  max resid 0.1245773 
Run 64 stress 0.1815262 
Run 65 stress 0.181948 
Run 66 stress 0.1917236 
Run 67 stress 0.1868432 
Run 68 stress 0.1828928 
Run 69 stress 0.1808185 
... New best solution
... Procrustes: rmse 0.006158855  max resid 0.03531337 
Run 70 stress 0.1819127 
Run 71 stress 0.1815252 
Run 72 stress 0.1886606 
Run 73 stress 0.186167 
Run 74 stress 0.180948 
... Procrustes: rmse 0.009886941  max resid 0.05098863 
Run 75 stress 0.1816035 
Run 76 stress 0.1935237 
Run 77 stress 0.1812551 
... Procrustes: rmse 0.01336276  max resid 0.08156214 
Run 78 stress 0.1822264 
Run 79 stress 0.1814582 
Run 80 stress 0.1814144 
Run 81 stress 0.1828653 
Run 82 stress 0.1823118 
Run 83 stress 0.1978016 
Run 84 stress 0.1810713 
... Procrustes: rmse 0.0369974  max resid 0.1455491 
Run 85 stress 0.1808181 
... New best solution
... Procrustes: rmse 0.0002833112  max resid 0.0008111158 
... Similar to previous best
*** Solution reached

T1_rare<-ggordiplots::gg_ordiplot(ord = ord_t1_rare, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T1")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)

T1_rare_plot<-T1_rare$plot + theme_classic() + xlim(-0.4, 0.4) + ylim(-0.4, 0.4)  + guides(color=guide_legend("Treatment")) + xlab("NMDS 1") + ylab("NMDS 2")
T1_rare_plot

library(cowplot)
my_legend <- get_legend(T1_rare_plot)
library(ggpubr)
as_ggplot(my_legend)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordinationlegend.pdf")

Saving 7.29 x 4.51 in image

T1_rare_plot<-T1_rare_plot + theme(legend.position = "none")
T1_rare_plot
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_rare_T1.pdf")

Saving 7.29 x 4.51 in image

ord_t2_rare<-ordinate(physeq = subset_samples(ps_rare, Time=="T2"), method = "NMDS", distance = "bray", k=3, trymax= 100)

Square root transformation
Wisconsin double standardization
Run 0 stress 0.1534021 
Run 1 stress 0.1495833 
... New best solution
... Procrustes: rmse 0.06591139  max resid 0.2432552 
Run 2 stress 0.1529359 
Run 3 stress 0.149299 
... New best solution
... Procrustes: rmse 0.01311101  max resid 0.05357392 
Run 4 stress 0.1498585 
Run 5 stress 0.153965 
Run 6 stress 0.1535223 
Run 7 stress 0.1501546 
Run 8 stress 0.1508771 
Run 9 stress 0.1500086 
Run 10 stress 0.154203 
Run 11 stress 0.1496264 
... Procrustes: rmse 0.04805853  max resid 0.2154931 
Run 12 stress 0.1546731 
Run 13 stress 0.1496216 
... Procrustes: rmse 0.04712721  max resid 0.2093185 
Run 14 stress 0.15183 
Run 15 stress 0.1496336 
... Procrustes: rmse 0.03080855  max resid 0.1376397 
Run 16 stress 0.1496321 
... Procrustes: rmse 0.03047098  max resid 0.13453 
Run 17 stress 0.1500846 
Run 18 stress 0.1495175 
... Procrustes: rmse 0.02694526  max resid 0.1366575 
Run 19 stress 0.1517768 
Run 20 stress 0.1498159 
Run 21 stress 0.1525828 
Run 22 stress 0.1496191 
... Procrustes: rmse 0.04687338  max resid 0.205771 
Run 23 stress 0.1558324 
Run 24 stress 0.1558559 
Run 25 stress 0.1517661 
Run 26 stress 0.1493063 
... Procrustes: rmse 0.00273002  max resid 0.01544446 
Run 27 stress 0.1543836 
Run 28 stress 0.1493533 
... Procrustes: rmse 0.02753152  max resid 0.1459909 
Run 29 stress 0.1497401 
... Procrustes: rmse 0.02124904  max resid 0.09580409 
Run 30 stress 0.1536117 
Run 31 stress 0.152943 
Run 32 stress 0.153367 
Run 33 stress 0.1502911 
Run 34 stress 0.1492996 
... Procrustes: rmse 0.009094635  max resid 0.05431342 
Run 35 stress 0.1500939 
Run 36 stress 0.157447 
Run 37 stress 0.1501923 
Run 38 stress 0.1503452 
Run 39 stress 0.1499623 
Run 40 stress 0.1565076 
Run 41 stress 0.1552341 
Run 42 stress 0.1493037 
... Procrustes: rmse 0.01222299  max resid 0.07016078 
Run 43 stress 0.1542397 
Run 44 stress 0.1504335 
Run 45 stress 0.1502563 
Run 46 stress 0.1493242 
... Procrustes: rmse 0.006378492  max resid 0.03651169 
Run 47 stress 0.1493052 
... Procrustes: rmse 0.002363587  max resid 0.01335982 
Run 48 stress 0.1496238 
... Procrustes: rmse 0.04801997  max resid 0.2142248 
Run 49 stress 0.1493332 
... Procrustes: rmse 0.007566515  max resid 0.04258916 
Run 50 stress 0.1508837 
Run 51 stress 0.1522037 
Run 52 stress 0.1492999 
... Procrustes: rmse 0.01056149  max resid 0.06267914 
Run 53 stress 0.1504227 
Run 54 stress 0.1495976 
... Procrustes: rmse 0.01568511  max resid 0.07203582 
Run 55 stress 0.154981 
Run 56 stress 0.1541634 
Run 57 stress 0.1492997 
... Procrustes: rmse 0.001503889  max resid 0.006022071 
... Similar to previous best
*** Solution reached

T2_rare<-ggordiplots::gg_ordiplot(ord = ord_t2_rare, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T2")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)

T2_rare_plot<-T2_rare$plot + theme_classic() + xlim(-0.4, 0.4) + ylim(-0.4, 0.4)+ theme(legend.position = "none")  + xlab("NMDS 1") + ylab("NMDS 2")
T2_rare_plot

Warning: Removed 3 rows containing missing values (geom_point).

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_rare_T2.pdf")

Saving 7.29 x 4.51 in image
Warning: Removed 3 rows containing missing values (geom_point).

#G166SG identified as outlier based on plots with it included. Removed to create plot. 
ps_rare <-  subset_samples(ps_rare, sample_names(ps_rare) != "G166SG")
ord_t3_rare<-ordinate(physeq = subset_samples(ps_rare, Time=="T3"), method = "NMDS", distance = "bray", k=3, trymax= 100)

Square root transformation
Wisconsin double standardization
Run 0 stress 0.1676017 
Run 1 stress 0.1674408 
... New best solution
... Procrustes: rmse 0.04160458  max resid 0.2358723 
Run 2 stress 0.1676017 
... Procrustes: rmse 0.04176362  max resid 0.2370166 
Run 3 stress 0.1676781 
... Procrustes: rmse 0.03928737  max resid 0.234611 
Run 4 stress 0.1674409 
... Procrustes: rmse 0.0002702533  max resid 0.00115854 
... Similar to previous best
Run 5 stress 0.1674412 
... Procrustes: rmse 0.0009690672  max resid 0.005153489 
... Similar to previous best
Run 6 stress 0.1674405 
... New best solution
... Procrustes: rmse 0.0005420413  max resid 0.002292957 
... Similar to previous best
Run 7 stress 0.1674412 
... Procrustes: rmse 0.00046538  max resid 0.002452812 
... Similar to previous best
Run 8 stress 0.1752381 
Run 9 stress 0.1674408 
... Procrustes: rmse 0.0002149259  max resid 0.0009482109 
... Similar to previous best
Run 10 stress 0.1674412 
... Procrustes: rmse 0.0005201256  max resid 0.002801695 
... Similar to previous best
Run 11 stress 0.1674417 
... Procrustes: rmse 0.0006125434  max resid 0.002951462 
... Similar to previous best
Run 12 stress 0.167681 
... Procrustes: rmse 0.03993784  max resid 0.2358893 
Run 13 stress 0.1674407 
... Procrustes: rmse 0.0003510634  max resid 0.001544427 
... Similar to previous best
Run 14 stress 0.1676786 
... Procrustes: rmse 0.03985713  max resid 0.235582 
Run 15 stress 0.1674409 
... Procrustes: rmse 0.0005786015  max resid 0.003398217 
... Similar to previous best
Run 16 stress 0.1676795 
... Procrustes: rmse 0.03999965  max resid 0.2356988 
Run 17 stress 0.1674408 
... Procrustes: rmse 0.000393707  max resid 0.001767988 
... Similar to previous best
Run 18 stress 0.167441 
... Procrustes: rmse 0.0003750548  max resid 0.001801765 
... Similar to previous best
Run 19 stress 0.1674412 
... Procrustes: rmse 0.0005055888  max resid 0.002397934 
... Similar to previous best
Run 20 stress 0.1676016 
... Procrustes: rmse 0.0417145  max resid 0.2369253 
*** Solution reached

T3_rare<-ggordiplots::gg_ordiplot(ord = ord_t3_rare, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T3")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)

T3_rare_plot<-T3_rare$plot + theme_classic() + xlim(-0.4, 0.4) + ylim(-0.4, 0.4)+ theme(legend.position = "none")  + xlab("NMDS 1") + ylab("NMDS 2")
T3_rare_plot

Warning: Removed 1 rows containing missing values (geom_point).

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_rare_T3.pdf")

Saving 7.29 x 4.51 in image
Warning: Removed 1 rows containing missing values (geom_point).

library(ggpubr)
ggarrange(T1_rare_plot, T2_rare_plot, T3_rare_plot, ncol = 1)

Warning: Removed 3 rows containing missing values (geom_point).
Warning: Removed 1 rows containing missing values (geom_point).

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_combined_rare_within_group_beta.pdf", width = 10, height = 10)

CAP ordination plots rarefied - not used

t1_dist <- distance(subset_samples(ps_rare, Time=="T1"), method="bray") #get wUnifrac and save
t1_table<-as.matrix(dist(t1_dist)) #transform wUnifrac index
ord_t1_rare_cap <- capscale(t1_table ~ Herbicide, data.frame(sample_data(subset_samples(ps_rare, Time == "T1"))))
T1_rare<-ggordiplots::gg_ordiplot(ord = ord_t1_rare_cap, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T1")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T1_rare$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_rare_T1_cap.pdf")


t2_dist <- distance(subset_samples(ps_rare, Time=="T2"), method="bray") #get wUnifrac and save
t2_table<-as.matrix(dist(t2_dist)) #transform wUnifrac index
ord_t2_rare_cap <- capscale(t2_table ~ Herbicide, data.frame(sample_data(subset_samples(ps_rare, Time == "T2"))))
T2_rare<-ggordiplots::gg_ordiplot(ord = ord_t2_rare_cap, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T2")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T2_rare$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_rare_T2_cap.pdf")


#G166SG identified as outlier based on plots with it included. Removed to create plot. 
ps_rare <-  subset_samples(ps_rare, sample_names(ps_rare) != "G166SG")
t3_dist <- distance(subset_samples(ps_rare, Time=="T3"), method="bray") #get wUnifrac and save
t3_table<-as.matrix(dist(t3_dist)) #transform wUnifrac index
ord_t3_rare_cap <- capscale(t3_table ~ Herbicide, data.frame(sample_data(subset_samples(ps_rare, Time == "T3"))))
T3_rare<-ggordiplots::gg_ordiplot(ord = ord_t3_rare_cap, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T3")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T3_rare$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_rare_T3_cap.pdf")
ggplot_build(T3_rare$plot)$data

Pairwise adonis testing

ps_pairwiseadonis<-function(physeq){
  otu_tab<-data.frame(phyloseq::otu_table(physeq))
  md_tab<-data.frame(phyloseq::sample_data(physeq))
    if(taxa_are_rows(physeq)== T){
       physeq_dist<-parallelDist::parDist(as.matrix(t(otu_tab)), method = "bray")}
            else{physeq_dist<-parallelDist::parDist(as.matrix(otu_tab), method = "bray")}
pairwiseAdonis::pairwise.adonis(x = physeq_dist, factors = md_tab$Herbicide, p.adjust.m = "none", perm = 1000)
}

ps_t1<-subset_samples(ps_rare_sub, Time == "T1")
ps_t1<-prune_taxa(taxa_sums(ps_t1) > 2, ps_t1)

ps_t2<-subset_samples(ps_rare_sub, Time == "T2")
ps_t2<-prune_taxa(taxa_sums(ps_t2) > 2, ps_t2)

ps_t3<-subset_samples(ps_rare_sub, Time == "T3")
ps_t3<-prune_taxa(taxa_sums(ps_t3) > 2, ps_t3)


ps_pairwiseadonis(ps_t1)

ps_pairwiseadonis(ps_t2)

ps_pairwiseadonis(ps_t3)

Pairwise betadispr by treatment, time and mode

ps_betadispr<-function(physeq, groupingvar = "Groupingvar"){
  otu_tab<-data.frame(phyloseq::otu_table(physeq))
  md_tab<-data.frame(phyloseq::sample_data(physeq))
    if(taxa_are_rows(physeq)== T){
       physeq_dist<-parallelDist::parDist(as.matrix(t(otu_tab)), method = "bray")}
            else{physeq_dist<-parallelDist::parDist(as.matrix(otu_tab), method = "bray")}
                mod<-vegan::betadisper(physeq_dist, md_tab[,groupingvar])
        ## Perform test
                print(anova(mod))
        ## Permutation test for F
                pmod <- vegan::permutest(mod, permutations = 1000, pairwise = TRUE)
                print(pmod)
                print(boxplot(mod))
}

permute test of disperson

ps_betadispr(subset_samples(ps_rare_sub, Time == "T1"), groupingvar = "Mode")
ps_betadispr(subset_samples(ps_rare_sub, Time == "T2"), groupingvar = "Mode")
ps_betadispr(subset_samples(ps_rare_sub, Time == "T3"), groupingvar = "Mode")


ps_betadispr(subset_samples(ps_rare_sub, Mode == "Chemical"), groupingvar = "Time")
ps_betadispr(subset_samples(ps_rare_sub, Mode == "Non-Treated"), groupingvar = "Time")
ps_betadispr(subset_samples(ps_rare_sub, Mode == "Hand"), groupingvar = "Time")


ps_betadispr(subset_samples(ps_rare_sub, Time == "T1"), groupingvar = "Herbicide")
ps_betadispr(subset_samples(ps_rare_sub, Time == "T2"), groupingvar = "Herbicide")
ps_betadispr(subset_samples(ps_rare_sub, Time == "T3"), groupingvar = "Herbicide")


ps_betadispr(subset_samples(ps_rare_sub, Herbicide == "Glyphosate"), groupingvar = "Time")
ps_betadispr(subset_samples(ps_rare_sub, Herbicide == "Atrazine-Mesotrione"), groupingvar = "Time")
ps_betadispr(subset_samples(ps_rare_sub, Herbicide == "Dicamba"), groupingvar = "Time")
ps_betadispr(subset_samples(ps_rare_sub, Herbicide == "Handweeded"), groupingvar = "Time")
ps_betadispr(subset_samples(ps_rare_sub, Herbicide == "Non-Treated"), groupingvar = "Time")

ps_betadispr(ps_rare_sub, groupingvar = "Herbicide")
ps_betadispr(ps_rare_sub, groupingvar = "Mode")
ps_betadispr(ps_rare_sub, groupingvar = "Time")

box and whisker plots of pairwise distance within group distances

#remotes::install_github("antonioggsousa/micrUBIfuns")
library(micrUBIfuns)
T1_beta<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T1"), method = "bray", group = "Herbicide")
T1_beta_plot <- T1_beta$plot
T1_beta_plot <- T1_beta_plot + theme_classic()+ guides(color=guide_legend("Treatment")) + ylab("Bray-Curtis Dissimilarity") + xlab("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.5, 0.75)
T1_beta_plot

Warning: Removed 1 rows containing non-finite values (stat_boxplot).
Warning: Removed 1 rows containing missing values (geom_point).

my_legend <- get_legend(T1_beta_plot)

Warning: Removed 1 rows containing non-finite values (stat_boxplot).
Warning: Removed 1 rows containing missing values (geom_point).

as_ggplot(my_legend)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_beta_legend.pdf")

Saving 7.29 x 4.51 in image

T1_beta_plot<-T1_beta_plot+ theme(legend.position = "none") 
T1_beta_plot

Warning: Removed 1 rows containing non-finite values (stat_boxplot).
Warning: Removed 1 rows containing missing values (geom_point).

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T1_rare_withingroup_beta.pdf")

Saving 7.29 x 4.51 in image
Warning: Removed 1 rows containing non-finite values (stat_boxplot).
Warning: Removed 1 rows containing missing values (geom_point).

T1_beta_df<- T1_beta$data
T1_betamod<-aov(formula = beta_div_value ~ group ,data = T1_beta_df)
summary(T1_betamod)

             Df  Sum Sq  Mean Sq F value   Pr(>F)    
group         4 0.02391 0.005978     5.4 0.000334 ***
Residuals   282 0.31223 0.001107                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

TukeyHSD(x = T1_betamod, which = "group")

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = beta_div_value ~ group, data = T1_beta_df)

$group
                                        diff          lwr          upr     p adj
Dicamba-Atrazine-Mesotrione     -0.001507071 -0.018186144  0.015172002 0.9991606
Glyphosate-Atrazine-Mesotrione   0.001896970 -0.015523754  0.019317693 0.9982525
Handweeded-Atrazine-Mesotrione  -0.023576431 -0.041939486 -0.005213376 0.0044684
Non-Treated-Atrazine-Mesotrione -0.013193939 -0.029873012  0.003485134 0.1934686
Glyphosate-Dicamba               0.003404040 -0.013275033  0.020083113 0.9805780
Handweeded-Dicamba              -0.022069360 -0.039730381 -0.004408340 0.0061800
Non-Treated-Dicamba             -0.011686869 -0.027589741  0.004216003 0.2601887
Handweeded-Glyphosate           -0.025473401 -0.043836456 -0.007110346 0.0015994
Non-Treated-Glyphosate          -0.015090909 -0.031769982  0.001588164 0.0971493
Non-Treated-Handweeded           0.010382492 -0.007278529  0.028043512 0.4896459

T2_beta<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T2"), method = "bray", group = "Herbicide")
T2_beta_plot <- T2_beta$plot
T2_beta_plot <- T2_beta_plot+ theme_classic() + theme(legend.position = "none") + ylab("Bray-Curtis Dissimilarity") + xlab("") + ggtitle("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.5, 0.75)
T2_beta_plot
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T2_rare_withingroup_beta.pdf")

Saving 7.29 x 4.51 in image

T2_beta_df<- T2_beta$data
T2_betamod<-aov(formula = beta_div_value ~ group ,data = T2_beta_df)
summary(T2_betamod)

             Df Sum Sq  Mean Sq F value   Pr(>F)    
group         4 0.0328 0.008212   5.471 0.000303 ***
Residuals   260 0.3902 0.001501                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

TukeyHSD(x = T2_betamod, which = "group")

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = beta_div_value ~ group, data = T2_beta_df)

$group
                                        diff          lwr          upr     p adj
Dicamba-Atrazine-Mesotrione      0.004587879 -0.015705828  0.024881585 0.9716343
Glyphosate-Atrazine-Mesotrione  -0.011421549 -0.032812994  0.009969896 0.5850476
Handweeded-Atrazine-Mesotrione   0.007709091 -0.012584615  0.028002797 0.8348219
Non-Treated-Atrazine-Mesotrione -0.022036364 -0.042330070 -0.001742657 0.0257655
Glyphosate-Dicamba              -0.016009428 -0.037400872  0.005382017 0.2426725
Handweeded-Dicamba               0.003121212 -0.017172494  0.023414919 0.9933115
Non-Treated-Dicamba             -0.026624242 -0.046917949 -0.006330536 0.0034255
Handweeded-Glyphosate            0.019130640 -0.002260805  0.040522085 0.1039333
Non-Treated-Glyphosate          -0.010614815 -0.032006260  0.010776630 0.6518171
Non-Treated-Handweeded          -0.029745455 -0.050039161 -0.009451748 0.0007049

T3_beta<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T3"), method = "bray", group = "Herbicide") 
T3_beta$plot #+ scale_color_manual(values = c("#F8766D", "#A3A500",  "#00BF7D", "#00B0F6", "#E76BF3")) + 
T3_beta_plot <- T3_beta$plot
T3_beta_plot <- T3_beta_plot + theme_classic()+ theme(legend.position = "none") + ylab("Bray-Curtis Dissimilarity") + xlab("") + ggtitle("")
T3_beta_plot <-T3_beta_plot + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.5, 0.75)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T3_rare_withingroup_beta.pdf")

Saving 7.29 x 4.51 in image
Warning: Removed 1 rows containing non-finite values (stat_boxplot).
Warning: Removed 1 rows containing missing values (geom_point).

T3_beta_df<- T3_beta$data
T3_betamod<-aov(formula = beta_div_value ~ group ,data = T3_beta_df)
summary(T3_betamod)

             Df Sum Sq  Mean Sq F value   Pr(>F)    
group         4 0.0477 0.011924   10.39 7.93e-08 ***
Residuals   261 0.2994 0.001147                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

TukeyHSD(x = T3_betamod, which = "group")

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = beta_div_value ~ group, data = T3_beta_df)

$group
                                         diff          lwr          upr     p adj
Dicamba-Atrazine-Mesotrione      0.0183969697  0.001410694  0.035383246 0.0263568
Glyphosate-Atrazine-Mesotrione  -0.0007666667 -0.018752976  0.017219643 0.9999574
Handweeded-Atrazine-Mesotrione   0.0276454545  0.010659179  0.044631731 0.0001130
Non-Treated-Atrazine-Mesotrione  0.0315000000  0.013513690  0.049486310 0.0000250
Glyphosate-Dicamba              -0.0191636364 -0.037864911 -0.000462362 0.0415656
Handweeded-Dicamba               0.0092484848 -0.008493102  0.026990072 0.6075977
Non-Treated-Dicamba              0.0131030303 -0.005598244  0.031804305 0.3068486
Handweeded-Glyphosate            0.0284121212  0.009710847  0.047113396 0.0003917
Non-Treated-Glyphosate           0.0322666667  0.012652605  0.051880729 0.0000917
Non-Treated-Handweeded           0.0038545455 -0.014846729  0.022555820 0.9798083

library(ggpubr)
ggarrange(T1_beta_plot, T2_beta_plot, T3_beta_plot, ncol = 1)

Warning: Removed 1 rows containing non-finite values (stat_boxplot).
Warning: Removed 1 rows containing missing values (geom_point).
Warning: Removed 1 rows containing non-finite values (stat_boxplot).
Warning: Removed 1 rows containing missing values (geom_point).

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_combined_rare_within_group_beta.pdf", width = 5, height = 10)

Examination of dissimliarity across time points by treatment and then again by all chemical treatments combined.

T1_beta_df$Time<-"T1"
T2_beta_df$Time<-"T2"
T3_beta_df$Time<-"T3"


beta_div_T1_T2_T3 <- rbind(T1_beta_df, T2_beta_df, T3_beta_df)

beta_anova<-function(data, Herbicide = "Herbicide"){
  df_sub<- data %>% filter(group == Herbicide)
  mod<-aov(beta_div_value ~ Time, data = df_sub)
  print(summary(mod))
  print(TukeyHSD(mod, "Time"))
  boxplot(df_sub$beta_div_value ~ df_sub$Time)
}

beta_anova(beta_div_T1_T2_T3, Herbicide = "Non-Treated")

             Df  Sum Sq Mean Sq F value   Pr(>F)    
Time          2 0.02547 0.01274   17.94 9.08e-08 ***
Residuals   163 0.11573 0.00071                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = beta_div_value ~ Time, data = df_sub)

$Time
             diff          lwr          upr     p adj
T2-T1 -0.01621212 -0.027718988 -0.004705255 0.0030431
T3-T1  0.01578788  0.003603568  0.027972190 0.0071575
T3-T2  0.03200000  0.019331357  0.044668643 0.0000000

beta_anova(beta_div_T1_T2_T3, Herbicide = "Handweeded")

             Df  Sum Sq  Mean Sq F value Pr(>F)  
Time          2 0.01713 0.008567   4.195 0.0168 *
Residuals   152 0.31041 0.002042                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = beta_div_value ~ Time, data = df_sub)

$Time
             diff           lwr        upr     p adj
T2-T1  0.02391582  0.0024158024 0.04541585 0.0251783
T3-T1  0.02231582  0.0008158024 0.04381585 0.0399463
T3-T2 -0.00160000 -0.0219967123 0.01879671 0.9811772

beta_anova(beta_div_T1_T2_T3, Herbicide = "Dicamba")

             Df  Sum Sq  Mean Sq F value Pr(>F)
Time          2 0.00273 0.001366   0.977  0.378
Residuals   173 0.24171 0.001397               
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = beta_div_value ~ Time, data = df_sub)

$Time
              diff         lwr         upr     p adj
T2-T1 -0.001274747 -0.01740797 0.014858477 0.9809505
T3-T1 -0.009002020 -0.02513524 0.007131204 0.3864971
T3-T2 -0.007727273 -0.02457788 0.009123330 0.5252045

beta_anova(beta_div_T1_T2_T3, Herbicide = "Atrazine-Mesotrione")

             Df  Sum Sq  Mean Sq F value   Pr(>F)    
Time          2 0.02773 0.013867   11.41 2.22e-05 ***
Residuals   173 0.21026 0.001215                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = beta_div_value ~ Time, data = df_sub)

$Time
              diff         lwr          upr     p adj
T2-T1 -0.007369697 -0.02308574  0.008346347 0.5100816
T3-T1 -0.028906061 -0.04395303 -0.013859094 0.0000310
T3-T2 -0.021536364 -0.03658333 -0.006489397 0.0025367

beta_anova(beta_div_T1_T2_T3, Herbicide = "Glyphosate")

             Df  Sum Sq  Mean Sq F value   Pr(>F)    
Time          2 0.02597 0.012985    14.9 1.33e-06 ***
Residuals   142 0.12374 0.000871                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = beta_div_value ~ Time, data = df_sub)

$Time
             diff         lwr          upr     p adj
T2-T1 -0.02068822 -0.03474249 -0.006633939 0.0018742
T3-T1 -0.03156970 -0.04562397 -0.017515421 0.0000012
T3-T2 -0.01088148 -0.02562173  0.003858768 0.1909546

#regroup all chemical treatments together and rerun betadiv calcs within group. 
sample_data(ps_rare)$Mode<-sample_data(ps_rare)$Herbicide

index <- c("Dicamba", "Glyphosate", "Atrazine-Mesotrione", "Handweeded", "Non-Treated")
values <- c("Chemical", "Chemical", "Chemical", "Handweeded", "Non-Treated")

sample_data(ps_rare)$Mode<- as.factor(values[match(sample_data(ps_rare)$Mode, index)])

#+ scale_color_manual(values = c("#FFA500", "#00B0F6", "#E76BF3")) 


T1_beta_chemical_combined<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T1"), method = "bray", group = "Mode")
T1_beta_chemical_combined_plot <- T1_beta_chemical_combined$plot 
T1_beta_chemical_combined_plot<- T1_beta_chemical_combined_plot + theme_classic() + guides(color=guide_legend("Treatment")) + ylab("Bray-Curtis Dissimilarity") + xlab("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.5, 0.75) + scale_color_manual(values = c("#FFA500", "#00B0F6", "#E76BF3")) 
T1_beta_chemical_combined_plot

Warning: Removed 2 rows containing non-finite values (stat_boxplot).
Warning: Removed 2 rows containing missing values (geom_point).

my_legend <- get_legend(T1_beta_chemical_combined_plot)

Warning: Removed 2 rows containing non-finite values (stat_boxplot).
Warning: Removed 2 rows containing missing values (geom_point).

as_ggplot(my_legend)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_beta_combined_legend.pdf")

Saving 7.29 x 4.51 in image

T1_beta_chemical_combined_plot<-T1_beta_chemical_combined_plot+ theme(legend.position = "none")
T1_beta_chemical_combined_plot

Warning: Removed 2 rows containing non-finite values (stat_boxplot).
Warning: Removed 2 rows containing missing values (geom_point).

T2_beta_chemical_combined<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T2"), method = "bray", group = "Mode")
T2_beta_chemical_combined_plot <- T2_beta_chemical_combined$plot 
T2_beta_chemical_combined_plot<- T2_beta_chemical_combined_plot + theme_classic()+ theme(legend.position = "none") + ylab("Bray-Curtis Dissimilarity") + xlab("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.5, 0.75) + scale_color_manual(values = c("#FFA500", "#00B0F6", "#E76BF3")) 
T2_beta_chemical_combined_plot

T3_beta_chemical_combined<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T3"), method = "bray", group = "Mode")
T3_beta_chemical_combined_plot <- T3_beta_chemical_combined$plot 
T3_beta_chemical_combined_plot<- T3_beta_chemical_combined_plot + theme_classic()+ theme(legend.position = "none") + ylab("Bray-Curtis Dissimilarity") + xlab("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.5, 0.75) + scale_color_manual(values = c("#FFA500", "#00B0F6", "#E76BF3")) 
T3_beta_chemical_combined_plot

Warning: Removed 1 rows containing non-finite values (stat_boxplot).
Warning: Removed 1 rows containing missing values (geom_point).

ggarrange(T1_beta_chemical_combined_plot, T2_beta_chemical_combined_plot, T3_beta_chemical_combined_plot, ncol = 1)

Warning: Removed 2 rows containing non-finite values (stat_boxplot).
Warning: Removed 2 rows containing missing values (geom_point).
Warning: Removed 1 rows containing non-finite values (stat_boxplot).
Warning: Removed 1 rows containing missing values (geom_point).

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_combined_rare_within_group_beta_chemical_combined.pdf", width = 5, height = 10)

T1_beta_df_chemical_combined <- T1_beta_chemical_combined$data
T2_beta_df_chemical_combined<- T2_beta_chemical_combined$data
T3_beta_df_chemical_combined<- T3_beta_chemical_combined$data

T1_beta_df_chemical_combined$Time<-"T1"
T2_beta_df_chemical_combined$Time<-"T2"
T3_beta_df_chemical_combined$Time<-"T3"

m1<-aov(beta_div_value ~ group, data = T1_beta_df_chemical_combined)
summary(m1)
TukeyHSD(m1, "group")
boxplot(beta_div_value ~ group, data = T1_beta_df_chemical_combined)


m2<-aov(beta_div_value ~ group, data = T2_beta_df_chemical_combined)
summary(m2)
TukeyHSD(m2, "group")
boxplot(beta_div_value ~ group, data = T2_beta_df_chemical_combined)

m3<-aov(beta_div_value ~ group, data = T3_beta_df_chemical_combined)
summary(m3)
TukeyHSD(m3, "group")
boxplot(beta_div_value ~ group, data = T3_beta_df_chemical_combined)


beta_div_chemical_combined_T1_T2_T3 <- rbind(T1_beta_df_chemical_combined, T2_beta_df_chemical_combined, T3_beta_df_chemical_combined)

beta_anova(beta_div_chemical_combined_T1_T2_T3, Herbicide = "Chemical")
beta_anova(beta_div_chemical_combined_T1_T2_T3, Herbicide = "Hand")
beta_anova(beta_div_chemical_combined_T1_T2_T3, Herbicide = "Non-Treated")

treatment to control

plotDistances = function(p, m, s, d) {

  # calc distances
  wu = phyloseq::distance(p, m)
  wu.m = melt(as.matrix(wu))
  
  # remove self-comparisons
  wu.m = wu.m %>%
    filter(as.character(Var1) != as.character(Var2)) %>%
    mutate_if(is.factor,as.character)
  
  # get sample data (S4 error OK and expected)
  sd = data.frame(sample_data(p)) %>%
    select(s, d) %>%
    mutate_if(is.factor,as.character)
  sd$Herbicide <- factor(sd$Herbicide, levels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))
  
  # combined distances with sample data
  colnames(sd) = c("Var1", "Type1")
  wu.sd = left_join(wu.m, sd, by = "Var1")
  
  colnames(sd) = c("Var2", "Type2")
  wu.sd = left_join(wu.sd, sd, by = "Var2")
  
  #remove this line to plot all comparisons. 
  #wu.sd = wu.sd %>% filter(Type1 == "Hand" | Type1 == "Non-Treated")
  
  # plot
  ggplot(wu.sd, aes(x = Type2, y = value)) +
    theme_bw() +
    geom_point() +
    geom_boxplot(aes(color = ifelse(Type1 == Type2, "red", "black"))) +
    scale_color_identity() +
    facet_wrap(~ Type1, scales = "free_x") +
    theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) + 
    ggtitle(paste0("Distance Metric = ", m))
  
}

a<-plotDistances(p = subset_samples(physeq= ps_rare, Time=="T1"), m = "bray", s = "Barcode_ID_G", d = "Herbicide")
a <- a + ggtitle("Time 1 Bray-Curtis Dissimlarities")
#ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T1_rare_allgroup_beta.pdf")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T1_rare_allgroup_beta_multicomparison.pdf")
b<-plotDistances(p = subset_samples(physeq= ps_rare, Time=="T2"), m = "bray", s = "Barcode_ID_G", d = "Herbicide")
b <-b + ggtitle("Time 2 Bray-Curtis Dissimlarities")
#ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T2_rare_allgroup_beta.pdf")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T2_rare_allgroup_beta_multicomparison.pdf")
c<-plotDistances(p = subset_samples(physeq= ps_rare, Time=="T3"), m = "bray", s = "Barcode_ID_G", d = "Herbicide")
c<- c + ggtitle("Time 3 Bray-Curtis Dissimlarities")
#ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T3_rare_allgroup_beta.pdf")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T3_rare_allgroup_beta_multicomparison.pdf")

library(ggpubr)
ggarrange(a, b, c, ncol = 1)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_combined_rare_allgroup_beta_multi_comparison.pdf", width = 7, height = 10)

Taxon abundance bar plot

#create super long color vector
col_vector <- c("#000000", "#FFFF00", "#1CE6FF", "#FF34FF", "#FF4A46", "#008941", "#006FA6", "#A30059",
        "#FFDBE5", "#7A4900", "#0000A6", "#63FFAC", "#B79762", "#004D43", "#8FB0FF", "#997D87",
        "#5A0007", "#809693", "#FEFFE6", "#1B4400", "#4FC601", "#3B5DFF", "#4A3B53", "#FF2F80",
        "#61615A", "#BA0900", "#6B7900", "#00C2A0", "#FFAA92", "#FF90C9", "#B903AA", "#D16100",
        "#DDEFFF", "#000035", "#7B4F4B", "#A1C299", "#300018", "#0AA6D8", "#013349", "#00846F",
        "#372101", "#FFB500", "#C2FFED", "#A079BF", "#CC0744", "#C0B9B2", "#C2FF99", "#001E09",
        "#00489C", "#6F0062", "#0CBD66", "#EEC3FF", "#456D75", "#B77B68", "#7A87A1", "#788D66",
        "#885578", "#FAD09F", "#FF8A9A", "#D157A0", "#BEC459", "#456648", "#0086ED", "#886F4C",
        
        "#34362D", "#B4A8BD", "#00A6AA", "#452C2C", "#636375", "#A3C8C9", "#FF913F", "#938A81",
        "#575329", "#00FECF", "#B05B6F", "#8CD0FF", "#3B9700", "#04F757", "#C8A1A1", "#1E6E00",
        "#7900D7", "#A77500", "#6367A9", "#A05837", "#6B002C", "#772600", "#D790FF", "#9B9700",
        "#549E79", "#FFF69F", "#201625", "#72418F", "#BC23FF", "#99ADC0", "#3A2465", "#922329",
        "#5B4534", "#FDE8DC", "#404E55", "#0089A3", "#CB7E98", "#A4E804", "#324E72", "#6A3A4C",
        "#83AB58", "#001C1E", "#D1F7CE", "#004B28", "#C8D0F6", "#A3A489", "#806C66", "#222800",
        "#BF5650", "#E83000", "#66796D", "#DA007C", "#FF1A59", "#8ADBB4", "#1E0200", "#5B4E51",
        "#C895C5", "#320033", "#FF6832", "#66E1D3", "#CFCDAC", "#D0AC94", "#7ED379", "#012C58",
        
        "#7A7BFF", "#D68E01", "#353339", "#78AFA1", "#FEB2C6", "#75797C", "#837393", "#943A4D",
        "#B5F4FF", "#D2DCD5", "#9556BD", "#6A714A", "#001325", "#02525F", "#0AA3F7", "#E98176",
        "#DBD5DD", "#5EBCD1", "#3D4F44", "#7E6405", "#02684E", "#962B75", "#8D8546", "#9695C5",
        "#E773CE", "#D86A78", "#3E89BE", "#CA834E", "#518A87", "#5B113C", "#55813B", "#E704C4",
        "#00005F", "#A97399", "#4B8160", "#59738A", "#FF5DA7", "#F7C9BF", "#643127", "#513A01",
        "#6B94AA", "#51A058", "#A45B02", "#1D1702", "#E20027", "#E7AB63", "#4C6001", "#9C6966",
        "#64547B", "#97979E", "#006A66", "#391406", "#F4D749", "#0045D2", "#006C31", "#DDB6D0",
        "#7C6571", "#9FB2A4", "#00D891", "#15A08A", "#BC65E9", "#FFFFFE", "#C6DC99", "#203B3C",

        "#671190", "#6B3A64", "#F5E1FF", "#FFA0F2", "#CCAA35", "#374527", "#8BB400", "#797868",
        "#C6005A", "#3B000A", "#C86240", "#29607C", "#402334", "#7D5A44", "#CCB87C", "#B88183",
        "#AA5199", "#B5D6C3", "#A38469", "#9F94F0", "#A74571", "#B894A6", "#71BB8C", "#00B433",
        "#789EC9", "#6D80BA", "#953F00", "#5EFF03", "#E4FFFC", "#1BE177", "#BCB1E5", "#76912F",
        "#003109", "#0060CD", "#D20096", "#895563", "#29201D", "#5B3213", "#A76F42", "#89412E",
        "#1A3A2A", "#494B5A", "#A88C85", "#F4ABAA", "#A3F3AB", "#00C6C8", "#EA8B66", "#958A9F",
        "#BDC9D2", "#9FA064", "#BE4700", "#658188", "#83A485", "#453C23", "#47675D", "#3A3F00",
        "#061203", "#DFFB71", "#868E7E", "#98D058", "#6C8F7D", "#D7BFC2", "#3C3E6E", "#D83D66",

        "#2F5D9B", "#6C5E46", "#D25B88", "#5B656C", "#00B57F", "#545C46", "#866097", "#365D25",
        "#252F99", "#00CCFF", "#674E60", "#FC009C", "#92896B")

phylumGlommed <- tax_glom(ps_rare, "Phylum")

#t1
phylumGlommed_herb_t1 <- merge_samples(subset_samples(physeq= phylumGlommed, Time=="T1"), group = "Herbicide")
phylumGlommed_herb_t1 <- transform_sample_counts(phylumGlommed_herb_t1, function(OTU) OTU/sum(OTU))
sample_data(phylumGlommed_herb_t1)$Herbicide <- factor(sample_data(phylumGlommed_herb_t1)$Herbicide, levels = c(1, 2, 3, 4, 5), 
       labels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))

plot_bar(phylumGlommed_herb_t1, x = "Herbicide", fill = "Phylum")  + theme_classic() + ggtitle("Proportional Taxon Abundances Time 1") +
theme(legend.position="bottom") + guides(fill=guide_legend(nrow=6)) + geom_bar(stat="identity") + theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) + 
scale_fill_manual(values = col_vector)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_Taxon_barplot_t1.pdf")

#t2
phylumGlommed_herb_t2 <- merge_samples(subset_samples(physeq= phylumGlommed, Time=="T2"), group = "Herbicide")
phylumGlommed_herb_t2 <- transform_sample_counts(phylumGlommed_herb_t2, function(OTU) OTU/sum(OTU))
sample_data(phylumGlommed_herb_t2)$Herbicide <- factor(sample_data(phylumGlommed_herb_t2)$Herbicide, levels = c(1, 2, 3, 4, 5), 
       labels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))

plot_bar(phylumGlommed_herb_t2, x = "Herbicide", fill = "Phylum")  + theme_classic() + ggtitle("Proportional Taxon Abundances Time 1") +
theme(legend.position="bottom") + guides(fill=guide_legend(nrow=6)) + geom_bar(stat="identity") + theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) + 
scale_fill_manual(values = col_vector)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_Pt1/Figures/16S_Taxon_barplot_t2.pdf")

#t3
phylumGlommed_herb_t3 <- merge_samples(subset_samples(physeq= phylumGlommed, Time=="T3"), group = "Herbicide")
phylumGlommed_herb_t3 <- transform_sample_counts(phylumGlommed_herb_t3, function(OTU) OTU/sum(OTU))
sample_data(phylumGlommed_herb_t3)$Herbicide <- factor(sample_data(phylumGlommed_herb_t3)$Herbicide, levels = c(1, 2, 3, 4, 5), 
       labels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))

plot_bar(phylumGlommed_herb_t3, x = "Herbicide", fill = "Phylum")  + theme_classic() + ggtitle("Proportional Taxon Abundances Time 1") +
theme(legend.position="bottom") + guides(fill=guide_legend(nrow=6)) + geom_bar(stat="identity") + theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) + 
scale_fill_manual(values = col_vector)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_Pt1/Figures/16S_Taxon_barplot_t3.pdf")

Combined herbicide and time bar plot for exploration

sample_data(ps_rare)$herb_time<-paste(sample_data(ps_rare)$Herbicide, sample_data(ps_rare)$Time, sep = "_")
ps_rare_for_barplot <- prune_taxa(taxa_sums(ps_rare) > 50, ps_rare)
plot_bar(ps_rare_for_barplot, x= "herb_time", fill = "Family") + scale_fill_manual(values = col_vector) + geom_bar(stat="identity")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_BarPlot_Herbicide_Time.pdf", width = 20, height = 11)

Linear modeling of abundant taxa.

Tax_glom_Subset <- function (physeq, y = "taxLevel", nreturns = "Number of returns"){
   ps_1<- tax_glom(ps_rare_sub, taxrank = y )
    myTaxa <- names(sort(taxa_sums(ps_1), decreasing = TRUE)[1:nreturns])
       ps_1_sub <- prune_taxa(myTaxa, ps_1)
  return(ps_1_sub)
}


ps_rare_family_top25<-Tax_glom_Subset(physeq = ps_rare, nreturns = 25, y = "Family")

myTaxa <- names(sort(taxa_sums(ps_rare), decreasing = TRUE)[1:25])
ps_rare_asv_top25 <- prune_taxa(myTaxa, ps_rare)


#explore top 25 taxa with plot bar
plot_bar(ps_rare_family_top25, x= "herb_time", fill = "Family") + scale_fill_manual(values = col_vector) + geom_bar(stat="identity")
plot_bar(ps_rare_family_top25, x= "Time", fill = "Family") + scale_fill_manual(values = col_vector) + geom_bar(stat="identity")
plot_bar(ps_rare_family_top25, x= "Herbicide", fill = "Family") + scale_fill_manual(values = col_vector) + geom_bar(stat="identity")


#write function to wrangle data prior to anova

abund_aov_wrangle <- function (physeq, y = "Tax_Level"){
  tax<-tax_table(physeq)[,y]
   meta<-data.frame(sample_data(physeq))
  counts<-data.frame(otu_table(physeq))
  rownames(counts) <- tax[,1]
  counts<-data.frame(t(counts))
   counts$Time <- meta$Time 
   counts$Herbicide <- meta$Herbicide 
  counts$Herb_time <- meta$herb_time 
  return(counts)
}            

test<-abund_aov_wrangle(ps_rare_family_top25, y = "Family")



mod_abund<-function(count_tab, IV = "Groups to be tested") {
   for(j in 1:length(unique(count_tab[,"Herbicide"]))){
         data <- count_tab %>% filter(Herbicide == unique(count_tab$Herbicide)[j])
           #change this to the number of returns from the tax_glom_subset function
   for (i in 1:25) { 
            mod <- aov(unlist(data[i]) ~ matrix(data[,IV])) 
            #sanity check
            #print(c(j, i))
   if(summary(mod)[[1]][["Pr(>F)"]][1] <= 0.05) {
            #print(summary(mod))
     print(c(as.character(unique(count_tab[,"Herbicide"]))[j], names(data)[i]))
              boxplot(unlist(data[i]) ~ unlist(data[IV]), main =paste(names(data[i]), as.character(unique(count_tab[,"Herbicide"]))[j]), xlab= "Time", ylab="Abundance") 
           }
         }
      }
    }


mod_abund(test, IV = "Time")